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Abstract. We present the stochastic model of the galactic cosmic ray (GCR) particles 
transport in the heliosphere. Based on the solution of the Parker transport equation we 
developed models of the short-time variation of the GCR intensity, i.e. the Forbush decrease 
(Fd) and the 27-day variation of the GCR intensity. Parker transport equation being the 
Fokker-Planck type equation delineates non-stationary transport of charged particles in the 
turbulent medium. The presented approach of the numerical solution is grounded on solving of 
the set of equivalent stochastic differential equations (SDEs). We demonstrate the method of 
deriving from Parker transport equation the corresponding SDEs in the heliocentric spherical 
coordinate system for the backward approach. Features indicative the preeminence of the 
backward approach over the forward is stressed. We compare the outcomes of the stochastic 
model of the Fd and 27-day variation of the GCR intensity with our former models established 
by the finite difference method. Both models are in an agreement with the experimental data. 


1. Introduction 

First stochastic equation (Langevin’s) was linked with the Newton’s principle [I]. From the 
beginning of the XX century, the stochastic approach became more useful for describing physical 
random processes. The stochastic differential equations (SDEs) in conjunction with various 
Monte Carlo technics are broadly used in many fields, like physics, finance, biology, chemistry, 
engineering, or management science. The main statistical characteristic is the representation 
of the solution of the Fokker-Planck type equation as a probability density distributions. 
Employment of probabilistic description with Monte Carlo simulations allow to reduce the 
solution of the partial differential equation (PDE) describing the analyzed phenomena to the 
integration of SDEs. 

We apply stochastic methodology to model the galactic cosmic rays (GCR) transport in the 
heliosphere. Algorithms used in the particle transport simulations are mainly based on finite 
difference methods (e.g. m)- Galerkin methods are used for solving time-dependent high order 
PDEs (e.g. m)- The homotopy perturbation method [1] was recently developed for the numerical 
solution of various linear and nonlinear PDEs. Also exists approach to the PDE solution 
grounded on particle methods, characterized by low numerical dispersion [5]. 

During the propagation through the heliosphere, GCR particles are modulated by the solar 
wind and heliospheric magnetic field (HMF). Modulation of the GCR is a result of action of 
four primary processes: convection by the solar wind, diffusion on irregularities of HMF, particles 
drifts in the non-uniform magnetic field and adiabatic cooling. Transport of the GCR particles 


in the heliosphere can be described by the Parker transport equation [7]: 

^ = V • • V/) -{v, + U).Vf + ^{V- tJ)^, (1) 

where / = /(r, R, t) is an omnidirectional distribution function of three spatial coordinates 
r = (r,d,(p}, particles rigidity R and time t; U is solar wind velocity, Vd the drift velocity, and 
is the symmetric part of the diffusion tensor of the GCR particles. 

Employing the stochastic approach to solving the Parker transport equation is not the latest 
idea mm)- However, majority of models presented in the literature are used to determine the 
simulated spectra and compare it with experimental observations carried out by space probes 
as. Voyagers, AMS, BESS and PAMELA (e.g. [in]-[Il])- In this paper, we present models in 
which we not only reproduce the proton spectra but also simulate the short time GCR variations 
i.e. the Eorbush decrease (Ed) and 27-day variation. Additionally we compare the results of the 
stochastic modeling with our previous well grounded models developed by using finite difference 
method (EDM) of solution the Parker transport equation (e.g. |15]-[17]1. To perform the reliable 
comparison between two models we consider the same coefficients and parameters in the baseline 
Parker transport equation. Eurthermore, the parameters included attaining the GCR variations 
from the model are obtained by approximation of the experimental data. 

The aim of our paper is twofold. The first is to compose a consistent mathematical model of the 
GCR transport in the heliosphere by means of the SDEs. The second is to employ the created 
model to simulate the short-term variations being in agreement with the experimental data. 
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Figure 1. The sample pseudoparticles trajectories within the heliosphere. 


2. Stochastic approach 

In order to model the GCR transport in the heliosphere using stochastic methods equivalent 
SDEs must be obtained. In the first step, Parker transport equation (Eq. [T]) should be converted 
to a general form of Fokker-Planck equation (FPE). Depending on the direction of integration 
FPE, can be expressed in two forms |18j : 
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Figure 2. Trajectories of the pseudoparticles initialized with rigidity 10 GV from position 
r = lAU, 6 = 90°, if = 180° for A>0 and A<0 solar magnetic cycle. 




Figure 3. The histograms for the particles rigidity and exit time for the pseudoparticles 
initialized with rigidity 10 GV from position r = lAU, 9 = 90°, (p = 180° for A>0 and A<0 
solar magnetic cycle. 
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Figure 4. Latitude vs. longitude distribution of simulated pseudoparticles (protons) for A > 0 
and A < 0 solar magnetic cycle. 


SDEs equivalent to Eqs. [2] and [3] has a form (e.g. [T8]l: 


df = Ai ■ dt + Bij ■ dW, 


( 4 ) 


where r is the individual pseudoparticle trajectory in the phase space and dWi is the Wiener 
process, commonly written as dWi = y/di ■ dwi, with dwi being the randomly fluctuating term 
with Gaussian distribution. 

To solve Eq. [5] in both cases (backward and forward), at the onset we initiate pseudoparticle 
at some starting point in space and time and integrate alongside the pseudoparticle trajectory 
until it reaches the boundary. Choosing the forward or backward approach we need to bear in 
mind the problem that we want to solve. In the forward integration, pseudoparticles start from 
diverse boundary points, being for the GCR particles the entrance to the heliosphere. After 
that, their trajectories are traced up to the target position, e.g. 1 AU. Thus, a high number 
of pseudoparticles has to be initialized in order of obtaining a robust statistic because plenty 
of them do not reach the target position. The backward integration is much more effective in 
the case of the GCR propagation in the heliosphere. In the backward approach, the number of 
’useless’ particles is reduced. Pseudoparticles start from point of interest (e.g. 1 AU) and are 
traced backward in time until crossing the heliosphere boundary (in this paper this boundary 
is assumed at 100 AU, Eig. [1]). The value of the particle distribution function f(r,R) for the 
starting point can be found as an average of fiisiR) value for pseudoparticles characteristics at 
the entry positions, /(r, R) = -^ J2n=i fLis{R)i where fiisiR) is the cosmic ray local interstellar 
spectrum taken as in m for rigidity R of the particle at the exit/entrance point. 

The Parker equation (Eq. [T|) in the 3-D spherical coordinate system (r, 9, ip) as time-backward 
EPE diffusion equation has a form; 
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Figure 5. Simulated galactic protons rigidity spectra at Earth with respect to the LIS (solid 
line) at 100 AU for two cases: the A > 0 (dotted lines) and the A < 0 (dashed line) polarity 
epochs. 


In our model we are using a full 3D anisotropic diffusion tensor of the GCR particles Kij = 
K^p + containing the symmetric and antisymmetric parts presented in [20]. 

The drift velocity of the GCR particles is realized as: Vd^i = [H]. The equivalent to Eq. [5] 

set of SDEs with matrix Bij, {i,j = r, 9, ip) has a following form (the same can be found in [22jl: 
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The Euler—Maruyama scheme is used in order to integrate Eqs. [6] backward in time. To 
solve Eqs. [6] in spherical coordinates system we are using boundary conditions having form: 
(fi < 0 ^ ifi = ipi+27T, ifi > 2 it ^ ifi = ipi — 2TT, 9i < 0 ^ 9i = 9i + 27T and 9i > n ^ 9i = vr —|0j|. 
The reflecting boundary is considered as the inner radial boundary, ^ = 0 at r = 0.001 AU. An 
empty heliosphere constitutes an initial condition: fi{0.0lAU < r < W0AU,9,ip,R,0) = 0, as 
was shown in [12]. Performed tests proved that the 3000 pseudoparticles were enough number 
to simulate the short time variations of GCR particles with rigidity 10 GV. 

Fig. [2] shows the trajectories of simulated pseudoparticles (galactic protons) for the positive 
(A > 0) (upper panels) and for the negative (A < 0) polarity epoch (bottom panels), initially 
injected with rigidity 10 GV from position r = lAU, 6 = 90°, ip = 180°. The left panels in 
Fig. [2| present the radial position of pseudoparticles versus backward time, while the right panels 




































Figure 6. Changes of the expected amplitudes of the Fd of the GCR intensity at the Earth 
orbit, for the rigidity of 10 and 20 GV based on the solutions of the Parker transport equation by 
SDEs and EDM in comparison with the GCR intensity registered by Haleakala, Potchefstroom 
and Moscow neutron monitors during the Ed in 18 March - 4 April 2002. 


show the position in latitude versus time. Eig. [3] illustrates the binned exit rigidity of simulated 
pseudoparticles and corresponding binned propagation times for pseudoparticles. The exit time 
is a bit shorter for the A > 0 than for the A < 0 because of the large drift speed in the polar 
regions in A > 0. This is in agreement with [8]. Eig.[l]displays latitude vs. longitude distribution 
of simulated pseudoparticles initiated at Earth orbit with rigidity 10 GV for A > 0 and A < 0 
cycles. Eig. 0 ] shows the different character of heliospheric transport depending on various drift 
patterns in the A > 0 and A < 0 cycles. In the A > 0 cycle, pseudoparticles are transported 
mainly toward higher latitudes, while in the A < 0 period pseudoparticles are drifting outward 
mainly along the neutral sheet at low latitudes (Eig. [21 right panel). 

Eig. [5] demonstrates simulated rigidity spectra at the Earth for the A > 0 (dotted line), and for 
the A < 0 (dashed line) polarity epoch. The solid line represents the unmodulated spectrum 
(LIS) [19] at 100 AU. However, for the rigidity greater than 1 GV difference is subtle, but still 
visible, being in an agreement with Zhang results (compare with |8|, Eig. 3). 

3. Model of the Forbush decrease of the GCR intensity 

We present the model of the recurrent Fd taking place due to established corotating 
heliolongitudinal disturbances in the interplanetary space. Gorotating interaction regions (CIR) 
passing the Earth gradually diminishes the diffusion at the Earth orbit, causing larger scattering 
of the GCR particles, and in effect fewer GCR particles reach the Earth. We simulate this 
process by the gradual decrease and then the increase of the diffusion coefficient at the Earth 
orbit with respect the heliolongitudes. The diffusion coefficient of cosmic ray particles has 
a form: = Kq ■ K{r) ■ K{R,iy), where Kq = 10^^cm^/s, K{r) = 1 + 0.5 • (r/lAU) and 

K{R, u) = R^~'^. The exponent v pronounces the increase of the HMF turbulence in the vicinity 
of space where the Fd is created (e.g. mm), and is taken as: v = 1 + 0.3sin{^p — 90°) for 
r < 30AU and 90° < (p < 270°. We assume the existence of the two dimensional spiral Parker’s 
heliospheric magnetic field B [23| implemented through the angle ip = arctan{—B^/Br) in the 
3D anisotropic diffusion tensor Kij of GCR particles [20]. We compare the results obtained by 
the solution of the Parker transport equation by SDEs with our previous method of solution by 
EDM [23] assuming the same changes of all included parameters. 

The expected changes of the GCR intensity for the rigidity of 10 and 20 GV during the simulated 
Fd in comparison with the profiles of the daily GGR intensities recorded by the three neutron 
monitors with different cut off rigidities in 18 March - 4 April 2002 presents Fig. [6l One can 
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Figure 7. Changes of the expected amplitudes of the 27-day variation of the GCR intensity at 
the Earth orbit for the rigidity of 10 GV based on the solutions of the Parker transport equation 
by SDEs and EDM in comparison with the GCR intensity registered by Moscow neutron monitor 
for the 7 September- 3 October 2007. 


see that the proposed models are in a good coincidence with the experimental data and, as is 
expected the amplitude of the Ed decreases for higher rigidities. Moreover, the model of the Ed 
obtained based on the solution of the SDEs allows to reflect the stochastic character of the GCR 
particles distribution in the heliosphere and analyze the pseudoparticle trajectory thorough the 
3D heliosphere (Fig. [1]), which is not possible based on the solution of the Parker transport 
equation by the EDM (e.g. US!). This feature can allow a thorough analysis of the trajectory of 
the GCR particles when encountering occasional shock waves accompanying the sporadic Fd. 

4. Model of the 27-day variation of the GCR intensity 

The recurrence of the GCR intensity connected with the solar rotation is commonly called the 
27-day variation, although the durations can slightly differ. The 27-day variation of the GCR 
intensity is connected with the heliolongitudinal asymmetry of the electromagnetic conditions 
in the heliosphere. The recent minimum of solar activity between solar cycles No. 23 and 
24 was quite exceptional. Recurrent variations connected with corotating structures (~ 27 
days), at the end of 2007 and almost for the whole year 2008 were clearly established in 
all solar wind and interplanetary parameters. Consequently, the 27-day variations of cosmic 
ray intensity were clearly visible in a variety of cosmic ray counts of neutron monitors (e.g., 
[26l IT6] 1 and space probes (e.g., m)- We present the model of the 27-day variation of 
the GCR intensity considering an individual period of solar rotation starting at 2007.09.07. 
As it was stated by [2^- m the 27-day variation of the GCR intensity in the minimum 
epochs is preferentially related to the heliolongitudinal asymmetry of the solar wind velocity. 

In the model we apply approximation of the in situ measurements of the solar wind 
speed as source of the 27-day variation of the GCR intensity described by the formula: 
U = Uo{l — 0.31sin{ip + 6.10) -|- 0.06szn(2(/? -|- 0.82) — 0.10szn(3(/i — 1.04)), Uq = AOOkm/s. 

The diffusion coefficient of cosmic ray particles has a form: = Kq ■ K{r) ■ K{R), where 

Kq = lO^^cm^/s, K{r) = 1 -|- 0.5 . (r/lAU) and K{R) = R^-^. Fig. [7] compares the results of 
solution of the Parker transport equation by SDEs with solution by FDM for the GCR particles 
with rigidity R=10 GV. The model of the 27-day wave of the GCR intensity obtained by the by 
SDEs and FDM at the Earth orbit (1 AU, 6 = 90^) is in agreement with the data of Moscow 
NM (Fig. [7]). 





























5. Conclusion 

• We presented the model of the Fd and the 27-day variation of the GCR intensity obtained 
based on the stochastic approach to the solution of the Parker transport equation. The 
modeling results are in a good agreement with the neutron monitors data. 

• The SDEs were integrated backward in time in the spherical coordinates applying the full 
3D anisotropic diffusion tensor. 

• We showed the excellent agreement between the stochastic results and the finite difference 
method results presented in our previous papers. 

• The models obtained based on the solution of the SDEs allow to reflect the stochastic 
character of the GCR particles distribution in the heliosphere. Additionally this approach 
allows to trace the pseudoparticle trajectory through the 3D heliosphere which is not 
possible based on the solution of the Parker transport equation by finite difference methods. 
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